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ABSTRACT 

We explore self-similar dynamical processes in a spherical isothermal self-gravitational 
fluid with an emphasis on shocks and outline astrophysical applications of such shock 
solutions. The previous similarity shock solutions of Tsai & Hsu and of Shu et al. may 
be classified into two types: Class I solutions with downstream being free-fall collapses 
and Class II solutions with downstream being Larson-Penston (LP) type solutions. 
By the analyses of Lou & Shen and Shen & Lou, we further construct similarity shock 
solutions in the 'semi-complete space'. These general shock solutions can accommo- 
date and model dynamical processes of radial outflows (wind), inflows (accretion or 
contraction), subsonic oscillations, and free-fall core collapses all with shocks in var- 
ious settings such as star-forming molecular clouds, 'champagne flows' in Hii regions 
around luminous massive OB stars or surrounding quasars, dynamical connection be- 
tween the asymptotic giant branch phase to the proto-planetary nebula phase with 
a central hot white dwarf as well as accretion shocks around compact objects such 
as white dwarfs, neutron stars, and supermassive black holes. By a systematic explo- 
ration, we are able to construct families of infinitely many discrete Class I and Class 
II solutions matching asymptotically with a static outer envelope of singular isother- 
mal sphere; the shock solutions of Tsai & Hsu form special subsets. These similarity 
shocks travel at either subsonic or supersonic constant speeds. We also construct twin 
shocks as well as an 'isothermal shock' separating two fluid regions of two different 
yet constant temperatures. 

Key words: hydrodynamics — stars: AGB and post AGB — stars: formation — 
planetary nebulae: general — H n regions — quasars: general 



1 INTRODUCTION 

Self-gravitational inflows (contractions or accretions or col- 
lapses) and outflows (expansions or winds) of an isothermal 
fluid with spherical symmetry have been investigated over 
the past several decades in various astrophysical contexts 
(e.g., star formations, supernova explosions, formation and 
evolution of galaxy clusters etc) . When a fluid system is suf- 
ficiently away from initial and boundary conditions, it may 
evolve into self-similar behaviours (e.g., Sedov 1959; Landau 
& Lifshitz 1959). Under isothermal condition and spherical 
symmetry, Larson (1969a,b) and Penston (1969a,b) indepen- 
dently found self-similar solutions, now commonly referred 
to as Larson-Penston (LP) type solutions. Shu (1977) ob- 
tained another class of similarity solutions containing the 
so-called 'expansion wave collapse solution' (EWCS) that 
describes an expanding region of core collapse from a sur- 



rounding molecular cloud. In addition to these solutions, 
Hunter (1977) derived a class of discrete similarity solutions 
in the 'complete solution space' from t — > — oo to t — > +oo. 
These solutions may represent clouds which are accreted to 
a central mass point as an inward propagating compression 
wave driven by external gas pressure. Whitworth & Sum- 
mers (1985) expanded the solution space by introducing a 
bounded two-parameter continuum {zo, wo} to the known 
discrete solutions. The two parameters zo and wo corre- 
spond to «o and mo in the solutions of Shu (1977), which 
characterize the asymptotic solution forms as t — > — oo and 
t — > +oo respectively. Hunter (1986) noted that the addi- 
tional bands of solutions of Whitworth & Summers (1985) 
involve weak discontinuities across the sonic critical line 
(Boily & Lynden-Bell 1995) and may be unstable to per- 
turbations (e.g., Lazarus 1981). Generalizing these earlier 
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results, Lou & Shen (2004) recently obtained a new class 
of self-similar isothermal solutions, referred to as 'envelope 
expansion with core collapse' (EECC) solutions, which are 
characterized by an interior free-fall collapse towards the 
core with an exterior envelope expanding at a constant ra- 
dial speed at large radii. A broad class of EECC solutions 
without crossing the sonic critical line can be readily ob- 
tained; meanwhile, a special class of infinitely many discrete 
EECC solutions crossing the sonic critical line twice analyt- 
ically can also be constructed. Boily & Lynden-Bell (1995) 
and Murakami et al. (2004) studied self-similar solutions for 
spherical collapses with radiative effects included. The work 
of Fatuzzo et al. (2004) treated protostellar collapse in a 
polytropic gaseous sphere and found self-similar solutions 
with nonzero inward speed at large radii (see earlier empha- 
sis of Shen & Lou 2004 on this point); in their model consid- 
eration, the infall rates are determined by the inward speed 
Voo, the overdensity A, the polytropic index F for the back- 
ground equation of state, and the polytropic index 7 for the 
dynamic equation of state; they did not explore polytropic 
solutions crossing the sonic critical line. 

Based on prior results, Tsai & Hsu (1995) found two 
classes of isothermal self-similar shock solutions referred to 
as Class I and II solutions which connects a static singular 
isothermal sphere (SIS) envelope to a free-fall solution and 
a LP solution respectively (see § 3.1.1). Shu et al. (2002) ex- 
panded Tsai & Hsu's LP similarity shock solution to model 
the so-called 'champagne flows' driven by a shock into Hll 
regions surrounding luminous massive OB stars (see § 3.1.2). 
Further extending the work of Tsai & Hsu (1995) and Shu et 
al. (2002), Shen & Lou (2004) introduced Class I shock so- 
lutions for dynamical evolution of young stellar objects and 
Class II solutions to model 'champagne flows' of Hll regions 
surrounding OB stars. More generally, they pointed out that 
shock similarity solutions can be matched with various types 
of upstream solutions (see § 3.1.3). 

In this paper, we construct self-similar isothermal shock 
solutions by exploring diverse possibilities. In § 2, we sum- 
marize the basic isothermal fluid equations, the self-similar 
transformation and the isothermal shock conditions. In ref- 
erence to previous similarity shock solutions, we derive in § 3 
an infinite number of discrete shock solutions for Class I and 
Class II categories by exploring the a — v phase diagram and 
present new similarity shock solutions: twin shock solution 
and two-temperature flows separated by shocks. In § 4, we 
discuss astrophysical applications of these shock solutions. 
In 8 5, we summarize and discuss our results. 



where p is the gas mass density and u is the radial flow 
speed. The radial momentum equation reads 



du 
di 



GM 



(3) 



a 2 dp 
dr p dr 

where G = 6.67 x 10" 8 dyn cm 2 g 2 is the gravitational 
constant, a = (p/p) 1 ' 2 is the isothermal sound speed and p 
is the gas pressure. The Poisson equation is automatically 
satisfied by equations Q and ©. 

The dimensionless independent similarity variable is 
x — r/(at) and the similarity transformations are 



p(r,t) = 
u(r,t) = 



a(x) 
AnGt 2 
- av(x) 



M(r,t) = 



a-i 
~G 



m(x) , 



a 2 (j> , 



(4) 



where a(x), m(x) and v(x) are dimensionless reduced vari- 
ables for mass density p(r,t), enclosed mass M(r,t) and ra- 
dial flow speed u(r,t), respectively; $ is the gravitational 
potential such that —d$/dr = —GM(r,t)/r 2 and 4>( x ) is 
its reduced form. These reduced variables depend only on x 
(e.g., Shu 1977; Lou & Shen 2004). 

Transformation an d equation Q lead to 

, . dm „ dm 
m + (v - 



.dm 



d.r 



which gives m(x) in terms of x, v and a, namely 
m(x) = x 2 a(x — v) . 



(•>) 



(6) 



The physical requirement of a positive m(x) corresponds to 
the solution regime of x — v > 0. 

By transformation an( A solution (JSJ in equations (2) 
and ©, we obtain a pair of coupled nonlinear ordinary dif- 
ferential equations (ODEs): 



[(*- 
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^ a dx 



x 



(x - v) 



(x — v) 



(7) 



(8) 



[see eqs (11) and (12) of Shu (1977) and eqs (9) and (10) in 
Lou & Shen (2004)]. In the two nonlinear ODEs, the critical 
line is characterized by (a; — v) 2 — 1=0 and a = 2/x. 

The analytical asymptotic solutions of ODEs and 
JSJ are shown below (see Lou & Shen 2004 for details). 

(i) For x — » +00, the leading order terms are 



V 



A_ 



Ax , 



(9) 



2 FORMULATION OF SHOCK CONDITIONS 

We begin with the fluid equations in the Eulerian form in the 
spherical polar coordinates (r, 9, </)). By spherical symmetry, 
the mass conservation is given by 



dM_ + u dM_ _ Q 
dt dr 



dM 
dr 



- Anr 2 p , (1) 



where M(r, t) is the enclosed mass within radius r at time 
t. The equivalent form of continuity equation Q is simply 



dp 10.2s 

l£ + ->Yr {r PU) = ° 



dt 



(2) 



where V and A are two independent constants, referred to 
as velocity and mass parameters, respectively. The envelope 
expansion approaches a constant speed V at large radii, 
(ii) For x — > + , the leading terms are 

v _ > - [ ^!!l£ ) , a —> ( i^L ) , m m (10) 



I m \ V 2 
\2x 3 ) 



describing a central free-fall collapse with a constant dimen- 
sionless parameter mo for central mass accretion rate (Shu 
1977; Lou & Shen 2004; Shen & Lou 2004; Yu & Lou 2005). 

For the central LP- type solutions as x — > + , the leading 
order terms are 



2 

r 



B 



Bx 3 /3 



(11) 
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with a finite reduced central mass density B. 

(iii) Along the sonic critical line, the two eigensolutions 
are respectively 

-v{x) = (1 - x,) + \J- — l\(x — x») H , 

a(x) = — - — (—-l)(x-x*) + --- (12) 
for type 1 solutions, and 

— v(x) = (1 — x*) —(a; — £») + • • • , 

a;* 

2 2 

afx) = ;r(ai — x. t ) + ■ ■ ■ (13) 

x, x% 

for type 2 solutions, as defined by Shu (1977) and used in the 
subsequent work [e.g. Hunter (1977); Whitworth & Summers 
(1985); Lou & Shen (2004); Shen & Lou (2004)]. These ana- 
lytical asymptotic solutions provide necessary starting con- 
ditions for a fourth-order Runge-Kutta integration scheme 
(Press et al. 1986) to obtain global solutions numerically. 
We derive various classes of numerical solutions from the 
two coupled nonlinear ODEs (J7J and JHJ globally. 

For an isothermal gas, cooling effects by dust grains al- 
low efficient radiation of compressional heat generated dur- 
ing a collapse. The energy conservation involves radiative 
losses (e.g., Courant & Friedrichs 1976; Spitzer 1978). Un- 
der the isothermal approximation, we need to consider con- 
servations of mass and momentum across a shock, namely 

Pd(ud —u 8 ) = p u (u u - u s ) , (14) 

ad pd + PdUd(u d - u s ) = a 2 p u + p u u u (u u - u s ) , (15) 

where u, a and p are the radial gas flow speed, the isother- 
mal sound speed and the mass density, respectively, and sub- 
scripts d (u) denote the downstream (upstream) of a shock, 
u s = adXsd = a u x su = T s jt is the outward moving speed 
of the shock with r 3 being the shock radius. Using dimen- 
sionless reduced variables v(x), x, and a(x) to replace u(r, t) 
and p(r,t), we derive the isothermal shock jump conditions 
in terms of v(x), x, and ot(x), namely 

ctd(vd — Xsd)a,d = a u (v u — x su )a u , (16) 

a.d 2 ctd + a,d 2 a>dVd(vd — x 3 d) = a u 2 a u + a u 2 a u v u (v u — x su ) ■ 

(17) 

For an isothermal sound speed ratio r = ad/a u and rx s d = 
^su ; we derive 

Vd~ Xsd~T(v u — Xsu) = (TV d -V u )(v u ~Xsu)(vd-Xsd) , (18) 

a d /a u = (v u - x su )/[r(v d - x sd )] (19) 
(Shen & Lou 2004). For r = 1, eqns JTHJ and lO become 
(v u - Xs){v d - x s ) = 1 (20) 

Ctd/ctu = (v u - Xs)/(v d - X s ) , (21) 

where x s = x su = x s d is the reduced shock speed. 



Our main goal is to solve coupled nonlinear ODEs J7J 
and (jHJ and to apply the jump conditions across various 
shocks to match with reasonable boundary conditions for 
constructing global similarity shock solutions. Conceptually, 
this would help to conceive various possible scenarios involv- 
ing shocks. With proper adaptations, this would lead to in- 
teresting physical models for various astrophysical systems. 



3 SELF-SIMILAR SHOCK SOLUTIONS 

We can construct various similarity solutions from nonlin- 
ear ODEs and © with shocks satisfying shock jump 
conditions 1181 1 and 11911 to match with different analytical 
asymptotic solutions. In general, these similarity solutions 
can describe radial outflows (expansions or winds), inflows 
(contractions or accretions), oscillations, free-fall collapses 
with one or more shocks. We explore the range of the inde- 
pendent similarity variable x in the 'semi-complete space', 
i.e., from x — > + to x — > +oo. By the invariant property for 
the time reversal of the equations, it is easy to establish the 
correspondence between solutions in the 'complete solution 
space' (Hunter 1977) and solutions in the 'semi-complete 
solution space' (Lou & Shen 2004). 

In this section, we first review earlier results of Tsai & 
Hsu (1995), Shu et al. (2002) and Shen & Lou (2004). We 
then derive two classes of infinitely many discrete shock so- 
lutions including the twin shock solutions. Finally, we study 
shock solutions with the shock separating the flows of two 
different yet constant temperatures. 

3.1 Previous Shock Solution Results 

We now review the similarity solution structure in reference 
to earlier results of Tsai & Hsu (1995), Shu et al. (2002), 
and Shen & Lou (2004). The common feature of these shock 
solutions is that they have been constructed with r = 1 and 
contain only one shock. As x — > + , these shock solutions 
can match with inner free-fall collapses across the sonic crit- 
ical line or with inner LP solutions, while as x — > +oo, they 
match with the analytical asymptotic solution (i) by equa- 
tion These solutions can be broadly classified by their 
behaviours near x — > + into two types: Class I whose down- 
stream solutions are free-fall collapses and Class II whose 
downstream solutions are LP-type solutions. 

3.1.1 Shock Solutions of Tsai & Hsu (1995) 

Shu (1977) obtained the EWCS (the dotted curve in Figs.[TJ 
with the mass density parameter A —* 2 + and the velocity 
parameter V = in equation @, and the solution touches 
the sonic critical line at x» = 1 with mo = 0.975. Tsai & 
Hsu (1995) considered a self-similar shock travelling into a 
static SIS envelope with r = 1, x s d — x su = x s , v u = and 
a u = 2/x 2 , such that the shock jump conditions I20H and 
12 1 > appear in the form of 

v d = Xs — , a d = 2 . (22) 

Xs 

They showed two specific examples for the two classes of 
similarity shock solutions: Class I solution given by heavy 
solid curves in the two panels of Fig.^and Class II solution 
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Figure 1. The reduced radial speed — v(x) in linear scale (upper 
panel) and the reduced mass density a(x) in logarithmic scale 
(lower panel) versus x in linear scale; the dashed curves corre- 
spond to the sonic critical line, the dotted curves correspond to 
the EWCS of Shu (1977), the heavy solid curve and the heavy 
dash-dotted curve of Tsai & Hsu (1995) correspond to the shock 
solutions matching with the inner free-fall collapses and with in- 
ner LP-type solutions to a static SIS envelope respectively, the in- 
ner free-fall solution crosses the sonic critical line at x* = 0.0554 
and the B parameter of the inner LP solution is B = 7.90; the 
light dash-dotted curves are solutions with asymptotic breezes at 
large x found by Shu et al. (2002) with the different curves from 
top to bottom corresponding to the parameter B = 2/3, 1, 1.67 
and 4 respectively in both upper and lower panels. 



given by heavy dashed-dotted curves in the two panels of 
Fig. This Class I solution crosses the sonic critical line at 
x, — 0.0554 and the shock location is at x s = 1.26 for an 
outgoing shock travelling at a constant speed of 1.26 times 
the sound speed a. As this Class I solution connects the in- 
ner free-fall collapse through an expansion into a static SIS 
envelope by a shock, the solution has the same properties of 
the EWCS such that the outer gas envelope remains static 
with a r~ 2 density profile at large radii and the central gas 
falls towards the core at a constant rate mo = 0.105 with a 
r -3/2 {jgjjgjty profile and a r _1,/2 velocity profile. The central 
mass accretion rate M — moo 3 /G in this similarity shock 
solution is a factor of ~ 0.1 that predicted by the EWCS 
of Shu (1977). For the Class II solution, a central expan- 
sion with a finite central density (LP solution) matches to a 
static SIS envelope across a shock; it shows a higher shock 
speed and strength. The shock location is at x s — 1.34 for 
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Figure 2. The dimensionless negative reduced radial speed — v(x) 
(top panel) and reduced mass density a(x) (bottom panel) versus 
x. The dashed lines in both panels are the sonic critical line; the 
light solid curves are illustrating examples of Class I solutions 
with the downstream solution crossing the sonic critical line at 
x t = 0.23 and the shock location x s at 0.43, 1.03, 1.37, 1.43, 
respectively; the heavy solid curves are illustrating examples of 
Class II solutions with the reduced density of the central core B = 
0.1 and the shock location x s at 1.30, 2.00 and 2.80, respectively. 



a constant shock expansion speed of 1.34 times the sound 
speed a and a reduced core density B = 7.9. 

3.1.2 Shock Solutions of Shu et al. (2002) 

Shu et al. (2002) extended the Class II solution of Tsai & Hsu 
(1995) and used these solutions to model 'champagne flows' 
of Hll regions surrounding massive OB stars. By varying 
the B parameter with fixed V — 0, they obtained a class of 
outflow solutions and found both shock speed and strength 
decreasing with increasing B. In Fig. Q we plot different 
solutions (dash-dotted curves) with different reduced central 
density parameter B — 2/3, 1, 1.67, 4, 7.90, corresponding 
to the mass parameter A = 1.09, 1.28, 1.52, 1.832, 2.0 
respectively, and also corresponding to the shock location 
x s = 2.06, 1.95, 1.80, 1.53, 1.34, respectively. As B -> 0+ 
(e.g., B — 10 -6 ) numerically, the mass parameter A also 
approaches + , and the reduced speed v(x) converges to an 
invariant form with an invariant fastest and strongest shock 
solution where the shock is located at x s = 2.56 (see the 
heavy solid curve in Fig.^J. 

3.1.3 Shock Solutions of Shen & Lou (2004) 

By constructing a new family of Class I similarity shock so- 
lutions to match with various asymptotic flows, Shen & Lou 
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Figure 3. The dimensionless reduced speed v(x) (top panel) and 
a scaled reduced mass density R(x) where R(x) = x 2 a(x)/B (bot- 
tom panel) versus x. The dashed lines are the sonic critical line; 
the solid curves are the family of Class II solutions in the invari- 
ant form with the reduced density of the central core B — ► 0+ 
(viz., B = 10 -6 ); the heavy solid line is the so-called 'champagne 
breeze' solution of Shu et al. (2002). 



(2004) modelled (light solid lines in Fig. |5J the dynamical 
evolution of young stellar objects such as the Bok globule 
B335 system to account better for the observationally in- 
ferred density and velocity profiles as well as the central 
mass accretion rate. The downstream is part of CP 1 solu- 
tion for an envelope expansion with core collapse (EECC; 
Shen & Lou 2004); this EECC solution crosses twice the 
sonic critical line (viz., x — v — 1 and a = 2/x) at a;* = 0.23 
and x* — 1.65 analytically. By choosing different shock lo- 
cations x s — 0.43, 1.03, 1.37, 1.43 as in Fig. [5] we readily 
construct various upstream solutions to match with the an- 
alytical asymptotic solutions at x — > Too (see Table 1 of 
Shen & Lou 2004). Comparing with Shu et al. (2002), Shen 
& Lou (2004) significantly broadened the Class II solutions 
(Fig-EJ for possible 'champagne flows' of Hll regions around 
luminous OB stars by allowing for flows at large x. The up- 
stream solutions in Shen & Lou (2004) are not limited to 
the static SIS envelope or 'champagne breeze' (V — 0); by 
adjusting shock locations (or equivalently, shock speed), it is 
possible to match with either asymptotic outflow (expansion 
or wind) solutions (V > 0) or asymptotic inflow (contraction 
or accretion) solutions (V < 0). 

3.2 Similarity Shocks into a Static SIS Envelope 

In Tsai & Hsu (1995) and Shu et al. (2002), two classes of 
solutions were derived to match with the static SIS envelope 
with V = and A = 2 or A < 2. Here, we further explore 
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Figure 4. The phase diagram of a versus v at the meeting point 
if = 0.5. Each open circle symbol indicates an integration from 
X* with its value shown explicitly and each asterisk symbol in- 
dicates an integration from the shock location x s with its value 
also shown explicitly. The solid and dotted curves intersect; here 
we only show four intersections, from right to left, corresponding 
to four pairs of {x,(l), x 3 } as {0.0554, 1.26}, {5.00 X 10" 5 , 1.04}, 
{6.50 X 10" 8 , 0.980} and {4.75 X 10" 4 , 0.640}, respectively. The 
trend of each curve suggests that there may exist an infinite num- 
ber of discrete Class I shock solutions matched with a static SIS 
envelope with A = 2. 
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Figure 5. Negative reduced radial speed —v(x) versus x for the 
four shock solutions identified in the a — v phase diagram of Fig. HI 
The straight dashed line is the sonic critical line. Numerals 1, 3, 
4, 2 marked along the four solid curves indicate the number of 
stagnation points corresponding to the shock solutions obtained 
in Fig.[I]from right to left. Key parameters are shown in the plot. 
All these shock solutions diverge, approaching — oo as x — » + ; 
such diverging behaviours are displayed more explicitly in Fig. |^_ 
below with the a;— axis in the logarithmic scale. 
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Table 1. Class I and II similarity shock solutions matched with a static SIS envelope 



description 




m 










v u 




node 


= 0.0544 
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1.26 
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1.26 
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9.42 x 10" 
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Columns 1 through 8 provide properties of solutions (x* for Class I solutions and B 
for Class II solutions); reduced central mass accretion rate mo for Class I solutions; the 
shock location x s ; the downstream velocity the downstream density ay; the upstream 
velocity v u , the upstream density a u and the number of stagnation points. 
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Figure 6. Enlarged portions of Class I shock solution curves near 
x — > + of Fig. |21 emphasizing the diverging and oscillatory be- 
haviours of this family of shock solutions as x — > + . The x— axis 
is shown in the logarithmic scale. The dashed line on the right 
is the sonic critical line. The undulatory profiles of the curves 
marked by numerals 2, 3, 4 represent self-similar subsonic radial 
oscillations with two, three, four stagnation points, respectively. 
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Figure 7. The phase diagram of a versus v at a chosen meeting 
point xp = 0.5. Each open circle symbol indicates an integration 
from x — » 0"*" [e.g., x = 10 — 5 to start solution II 11 1 with the 
central reduced mass density B shown explicitly and each asterisk 
symbol indicates an integration from the shock location x s with 
its value marked explicitly. The solid and dotted curves intersect; 
the first four intersections, from right to left, correspond to the 
parameter pair {B, x 3 } as {7.90, 1.335}, {9.50xl0 4 , 1.04}, {l.lOx 
10 7 , 0.966} and {930, 0.633}, respectively. The variation trend of 
the phase curves suggest an infinite number of discrete Class II 
shock solutions matched with a static SIS envelope A = 2. 



Class I and Class II solutions matched with the static SIS en- 
velope by searching possible solutions in the phase diagram 
of v and a. This procedure, routinely tested earlier, was 
first introduced by Hunter (1977; see Lou & Shen 2004 for 
the EECC solutions). The main technical difference is that 
in the prior analyses the adjustable parameters are points 
along the sonic critical line [i.e., horizontal coordinates a;*(l) 
and x*(2) as in Lou & Shen (2004)], while in our cases, more 
parameters are involved including x*(l), the shock location 
x s , the mass parameter A and the velocity parameter V = 0. 
We take one parameter as fixed and adjust the other two to 
match solutions in the phase diagram of v and a. 

We first focus on Class I shock solutions with a static 
SIS envelope, downstream solutions diverging at the centre, 



and a mass parameter A — 2 for upstream solutions. Tsai & 
Hsu (1995) derived a solution in this case and we will explore 
more thoroughly the phase diagram of v and a to construct 
more similarity shock solutions. In Fig.Q] we display the rel- 
evant phase diagram of v and a where we choose a meeting 
point at xp = 0.5 and integrate from x*(l) forward with a 
type 2 eigensolution solution to if = 0.5. Meanwhile, we 
use the shock jump condition to determine the v,i and as 
the starting condition to integrate backward from a chosen 
x s to if = 0.5. For i > i s , we simply set v = and A — 2 
for a static SIS envelope. Based on Fig. 21 in addition to the 
solution of Tsai & Hsu, three more similarity shock solutions 
belonging to this class can be found and are shown for — v 
versus x in Fig.|S|(the enlarged versions for small x are shown 
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Figure 8. Negative reduced speed —v(x) versus x of the four 
shock solutions for the four intersections in the a — v phase dia- 
gram in Fig.[7J The straight dashed line is the sonic line. Numerals 
0, 2, 3, 1 along the four solid curves denote the number of stag- 
nation points in the four shock solutions determined in Fig. |jj 
from right to left. Downstream portions are LP- type solutions. 
Key parameters are shown in the plot. 



in logarithmic scale in Fig.^J. The variation trend of spiral 
path in the phase diagram suggests that as a class, there ex- 
ists an infinite number of discrete similarity shock solutions 
that connect inner free-fall collapses to an outer static SIS 
envelope. Each of such solutions is expected to be uniquely 
associated with a number of stagnation points. From the 
intersections of the dotted and solid curves in Fig. |1] we ob- 
tain four shock solutions, from right to left, corresponding 
to the displayed solutions marked by 1, 3, 4, 2 in Figs.|S|and 
||3 these numerals denote the number of stagnation points, 
and the relevant parameters {a;*(l), x s } are {0.0554, 1.26}, 
{5.00 x 10" 5 ,1.04}, {6.50 x 10~ 8 , 0.980} and {4.75 x 10~ 4 , 
0.640}, respectively. The solution containing only one stag- 
nation point is just the solution first found by Tsai & Hsu 
(1995). We find three new solutions by exploring the phase 
diagram of v and a. In principle, we can construct more 
solutions following the sequence with increasing number of 
stagnation points. Behaviours of these solutions at x — > + 
are displayed in Fig. with velocity profiles crossing the 
line v = at different x values. As the value of becomes 
smaller, the number of stagnation points increases accord- 
ingly. For example, the shock in the shock solution with 
two stagnation points is an accretion shock such that both 
downstream and upstream fluids move towards the centre. 

We now turn to the Class II shock solutions that con- 
nect downstream LP solutions with a static SIS envelope. 
Tsai & Hsu (1995) first obtained a solution in this case; 
we shall demonstrate by examples the existence of more 
such Class II solutions with increasing number of stagna- 
tion points. Fig. [7J is a relevant phase diagram of v and a 
to identify shock solution matches. Here, we choose again 
a meeting point at xf — 0.5 and integrate forward from 
x — > + (e.g., x = 10 -5 ) where the analytical asymptotic 
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Figure 9. The phase diagram of v and a at a chosen meeting 
point xp = 3 for the shock solution with x* = 0.5. Each asterisk 
symbol indicates an integration from x s with its value shown 
explicitly and each open circle denotes an integration from x — » 
+oo (actually from x = 20) with the mass density parameter 
A also marked explicitly. From the intersection of the solid and 
dashed curves, we immediately determine the shock location x s 
and the mass density parameter A for this shock solution. 
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Figure 10. The Class I similarity shock solutions with V = 
0. With one stagnation point, there is only one solution 
matched with a static SIS envelope. For different x*, from 
left to right, the solid curves represent solutions of x* = 
0.0544, 0.10, 0.23, 0.50, 0.80. Note that the shock locations x a do 
not vary regularly by shifting x*. Details of this class of solutions 
are summarized in Table |21 

LP solution 1111 is imposed with different values of the cen- 
tral reduced density B. Meanwhile, we use the shock jump 
condition to determine Vd and as the starting condition 
to integrate from x s back to if = 0.5. For x > x a , we simply 
set v — and ^4 = 2. The variation trend of spiral structure 
in Fig. [7] is qualitatively similar to that in Fig. 2] pointing 
to a class of infinitely many discrete similarity shock so- 
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Table 2. Class I shock solutions with an asymptotic breeze 





m 


A 


Xs 




Old 


v u 




0.0554 


0.105 


2.00 


1.26 


0.463 


2.00 





1.26 


0.10 


0.188 


1.92 


1.35 


0.572 


1.72 


0.102 


1.11 


0.23 


0.404 


1.87 


1.37 


0.529 


1.54 


0.181 


1.09 


0.50 


0.738 


1.92 


1.23 


0.315 


1.63 


0.145 


1.38 


0.80 


0.939 


1.99 


1.07 


0.0841 


1.87 


0.0496 


1.81 



Columns 1 through 8 list the values of ic*(l) where the down- 
stream solution crosses the sonic critical line, the reduced central 
mass accretion rate mo , the mass density parameter A for the up- 
stream solution, the shock location x s , the downstream velocity 
Vd, the downstream density a^, the upstream velocity v u and the 
upstream density a u , respectively. 



lutions that connect inner LP solutions to an outer static 
SIS envelope. Based on the phase diagram of Fig. |H| the 
four intersections of the dotted and solid curves give rise to 
four Class II shock solutions, from right to left, correspond- 
ing to solutions marked by the numerals 0, 2, 3, 1 in Fig. |H| 
these numerals represent uniquely the number of stagnation 
points to identify solutions with corresponding parameters 
{B,x a } given by {7.90, 1.34}, {9.50 xlO 4 , 1.04}, {1.10 xlO 7 , 
0.966} and {930, 0.633}, respectively; the solution with no 
stagnation point is just the solution found by Tsai & Hsu 
(1995). As x —> + , the self-similar oscillation behaviours 
of downstream solutions of these Class II shock solutions 
are similar to those of Class I shock solutions, and as the 
central reduced density B becomes larger, the number of 
stagnation points increases. The main difference is that as 
x — > + , Class II shock solutions remain finite while Class I 
shock solutions diverge for core collapses. 

We now summarize the main results of this section. 
By numerical exploration and analysis, we derived discrete 
shock solutions matched with a static SIS envelope of V — 
and A — 2 within the semi-complete space of < x < +oo 
(i.e., for x > a; s , we set v = and A = 2). Given this SIS en- 
velope as the upstream condition, we explore classes of infi- 
nite number of discrete solutions for both Class I and Class II 
solutions whose downstreams are free-fall collapses and LP 
solutions, respectively. For Class I solutions as x — > + , the 
radial velocity and density profiles have asymptotic diverg- 
ing behaviours (1101 . while for Class II solutions as x — » 0, 
the radial velocity and density profiles have finite asymptotic 
behaviours 1 1 U characterized by a central reduced density 
B. For comparison, the behaviours x — » + of this fam- 
ily of shock solutions with an outer static SIS envelope are 
similar to those of EECC solutions, viz., both have sub- 
sonic oscillations in a self-similar manner as x* — > 0+. For 
Class I solutions as the value of x* decreases, the number of 
stagnation points where v = along the x— axis increases, 
while for Class II solutions as the value of reduced density 
B increases, the number of stagnation points increases. The 
Class I and Class II shock solutions shown by Tsai & Hsu 
(1995) which contain only one and no stagnation point along 
x— axis, respectively, are special examples belonging to these 
two classes of shock solutions. 



3.3 Similarity Shock Solutions of Breeze 

In this section, we allow mass parameter A to vary (i.e., 
4/2) for a fixed x* corresponding to a specific reduced 
central mass accretion rate mo. Parameters A and x s are 
adjusted to match the upstream and downstream solutions 
in the phase diagram of v and a. We integrate from a fixed 
x, (e.g., 0.5 as shown in Fig- El with a type 2 eigensolution 
to the shock location x s and get the Vd and ay- Across the 
shock, we get v u and a„ as the starting condition to integrate 
forward to a chosen meeting point xy = 3; meanwhile, we 
integrate back from upstream x — * +oo by varying the mass 
parameter A with a fixed speed parameter V — 0. 

Here, we apply the analytical asymptotic solution (i) as 
given by equation © at x — 20 as the starting condition 
to integrate back to xf=3. In Fig. [§] we take x* = 0.5 as 
an example to illustrate the procedure of determining the 
shock location and the mass parameter A for the upstream 
solution in the phase diagram of v and a. By doing so, we 
derive a new family of Class I shock solutions. For exam- 
ple, we pick values of x, as 0.0544, 0.1, 0.23, 0.5 and 0.8, 
corresponding to values of m = 0.105, 0.188, 0.404, 0.738, 
and 0.939, respectively. We find that for Class I shock solu- 
tions the value of shock velocity and strength and the mass 
parameter A do not vary with x* regularly as do those for 
Class II shock solutions. Class I shock solutions with differ- 
ent x* are shown in terms of — v and x in Fig. E3 in Table 
[5] we summarize these solution results. For example, the 
downstream solutions with a;*=0.1 and 0.5 can match with 
the upstream solutions for the same mass parameter A. 

To sum up, we present in this section the similarity so- 
lutions with asymptotic breezes, for which, we derived a se- 
ries of Class I solutions with different central mass accretion 
rate mo corresponding to different values of x*. This class 
of solutions can be applied to processes of cloud collapse in 
star formation. One important feature of these solutions is 
that the value of mo is adjustable within the range of to 
0.975 corresponding to different central mass accretion rates 
or shock locations. In one aspect, this family of solutions is 
analogous to the Class II shock solutions obtained by Shu et 
al. (2002), i.e., both the upstream solutions are breeze solu- 
tions. Two major differences are: (1) our solutions are Class 
I shock solutions with central free-fall collapses (i.e., diverg- 
ing solutions as x — » + ), while the downstream solutions of 
Shu et al. (2002) are LP solutions. (2) The shock location x s 
and the mass parameter A for upstream solutions vary reg- 
ularly as the value of reduced central density B decreases. 
In contrast, there is no obvious trend of regularity for our 
Class I shock solutions as x* decreases. 



3.4 Self-Similar Twin Shock Solutions 

In all previous studies, the similarity shock solutions con- 
tain only one shock and the Class I solutions crosses the 
sonic critical line only once analytically. In this section, we 
construct similarity shock solutions with twin shocks. 

In § 3.2, we find a Class I shock solution crossing the 
sonic critical line at x„ = 4.75 x 10~ 4 with a shock loca- 
tion x s = 0.640 and two stagnation points (see Fig. |SJ, 
and a Class II shock solution with a central reduced den- 
sity B — 930 and a shock location at x s = 0.633 (Fig. |HJ 
containing only one stagnation point (n.b., the one of v — » 
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Table 3. Class I and II twin shock solutions for different a;* (2) 





Description 




Z s (2) 


A 


T/ 

V 


0.80 


= 3.40 x 10~ 4 


0.670 


0.85 


1.43 


-0.528 




m = 6.20 x 10~ 4 












B = 1300 


0.660 


0.90 


1.52 


-0.430 








1.00 


1.77 


0.192 








1.07 


1.99 











1.10 


2.12 


0.109 








1.15 


2.36 


0.294 


0.90 


= 4.40 x 10~ 4 


0.640 


0.95 


1.77 


-0.198 




m = 8.18 x 10~ 4 












B=1000 


0.627 


1.00 


1.91 


-0.0708 








1.05 


2.10 


0.0831 



Columns 1 through 6 list values of x* (2) where the Class I twin 
shock solutions cross the sonic critical line the second time, the 
key parameters [a:*(l) and mo for Class I and B for Class II so- 
lutions, respectively], the first shock location rr s (l), the second 
shock location x s (2), the mass density parameter A for the up- 
stream solution and the speed parameter V, respectively. 
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Figure 11. The phase diagram of a versus v at a chosen meeting 
point = 0.5. Each solid dot denotes an integration from a;*(l) 
(its value shown explicitly) with type 2 eigcnsolution, each cross 
symbol denotes an integration from x — * + (e.g., x = 10 -5 ) with 
the value of the central reduced mass density B shown explicitly, 
and each open circle and asterisk symbols denote integrations 
from x*(2) = 0.8 and a;* (2) = 0.9 respectively via shock condi- 
tion to the meeting point xp = 0.5; numerals indicate the shock 
location x s . From the curve intersections, we determine the Class 
I and II shock solutions for x t = 0.8 and 0.9, respectively. 



as x — > + is not counted here). The portions of upstream 
solutions between the shock and the static SIS envelope of 
both shock solutions are tangent to the sonic critical line 
at the point of x*(2) = 1. In fact, the upstream solution is 
part of the EWCS (Shu 1977). As we shift the x*(2) location 
leftward into the interval < x < 1, the upstream solution 
of the first (i.e., leftmost) shock may cross the sonic critical 
line analytically with a type 2 eigensolution. We can then 
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Figure 12. Class I shock solutions with x*(2) = 0.8 and 0.9 
for —v{x) versus x denoted by the light and heavy solid curves, 
respectively. The straight dashed line is the sonic critical line. The 
two corresponding values of (1) are 4.40 xlO" 4 and 3.40 xl0~ 4 . 
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Figure 13. Class II shock solutions with x t = 0.8 and 0.9 for 
—v(x) versus x denoted by the light and heavy solid curves, re- 
spectively. The straight dashed line is the sonic critical line. Key 
parameters are shown in the plot. 



construct Class I and Class II twin shock solutions in the 
a — v phase diagram given this real possibility. 

Let us first focus on Class I twin shock solutions. We fix 
the value of x*(2) and adjust the value of x*(l) and the shock 
location x s (l) between x*(l) and x,(2) to find the similarity 
shock solution crossing the sonic critical line twice analyt- 
ically based on intersections in the a — v phase diagram. 
On this ground, it is possible to further construct a Class 
I similarity solution with twin shocks and an appropriate 
asymptotic solution at x — ♦ +oo. We take a meeting point 
at xf — 0.5 and integrate an upstream solution of the first 
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Figure 14. Class I and II twin shock solutions with i, = 0.8 rep- 
resented by the solid and dotted curves, respectively. The straight 
dashed line is the sonic critical line. Once the upstream solution 
of the first shock at x s (l) crosses the sonic critical line analyti- 
cally at x* = 0.8 with a type 2 eigensolution, we can choose the 
second shock location x 3 (2) = 0.85, 0.90, 1.00, 1.07, 1.10, 1.15 to 
match with various asymptotic solutions at x — > +oo. 
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Figure 15. Class I and II twin shock solutions with x* = 0.9 rep- 
resented by the solid and dotted curves, respectively. The straight 
dashed line is the sonic critical line. Once the upstream solution 
of the first shock at x s (l) crosses the sonic critical line analyti- 
cally, we choose for example x g (2) = 0.95, 1.00, 1.05 for the second 
shock location to match with asymptotic solution at x — > +oo. 



shock from x* (2) for Class I solutions by using a type 2 eigen- 
solution as a starting condition to a chosen shock location 
x 3 (l). Based on the velocity and density upstream of this 
point {v u , a u }, we use the shock jump conditions 12(11 and 
12111 to determine the velocity and density {vd, a<i} down- 
stream of x s (l) and to integrate from x a (l) to if = 0.5; 



meanwhile, we integrate forward from using another 

type 2 eigensolution to if = 0.5. In Fig. 1111 we show the 
phase diagram for the two cases of x*(2) = 0.8 and 0.9, re- 
spectively. From the intersections of the solid curve and the 
dotted curve as well as the dashed curve, we obtain the rele- 
vant parameter pair {x,(l) x s (l)} as {3.40 xlO -4 , 0.670} for 
the solution with a;, (2) = 0.8 and as {4.40 x 10~ 4 , 0.640} for 
the solution with £E*(2) = 0.9. We show the corresponding 
solutions for — v versus x in Fig. 1121 

For Class II solutions, we match the upstream and 
downstream solutions in the a — v phase diagram using 
a search procedure similar to that of finding the Class I 
twin shock solutions. Given a fixed £*, we adjust the value 
of central reduced mass density B and the shock location 
x 8 (l) in order to find matches. At a chosen meeting point 
x-p = 0.5, we first obtain different parameter sets {v, a} 
denoted by cross symbols for different B values by integrat- 
ing forward from x — » + with the LP solution. We then 
integrate backward from x, to an adjustable shock location 
x s (l); by imposing the shock condition, we continue to in- 
tegrate backward towards xp = 0.5 to obtain a pair of {v, 
a}. In Fig. 1111 we thus obtain open circle symbols for the 
case of x, — 0.8 with different x 3 (l) and asterisk symbols 
for the case of x* = 0.9 with different x 3 (l). From the in- 
tersections of the dash-dotted curve with the dotted and 
dashed curves, we find the parameter pair {B, x a (l)} as 
{1300, 0.660} and {1000, 0.627} for the shock solutions of 
x* = 0.8 and x, = 0.9, respectively. 

By matching the downstream and upstream solutions in 
the phase diagram, we obtain the solution structure around 
the first shock x a (l). For a solution passing the sonic crit- 
ical line at i, the second time, we can further construct 
the second shock by choosing the different location x a (2) of 
the second shock. For x, — 0.8, we choose x 3 (2) — 0.85, 
0.90, 1.00, 1.07, 1.10 and 1.15 while for x* = 0.9, we choose 
x s (2) = 0.95, 1.00, 1.05. The upstream solutions across this 
second shock can match with the analytical asymptotic so- 
lution J2J at x — > +00; the property of these upstream solu- 
tions can evolve from the inflow (contraction or accretion) 
to outflow (expansion or wind or breeze) , as the value of the 
shock location x 3 (2) increases. 

The main properties of the similarity twin shock solu- 
tions are: (1) there exist two stagnation points for the Class I 
twin shock solutions, while there is only one stagnation point 
for the Class II twin shock solutions excluding x — > + ; the 
spherical stagnation surfaces of zero flow speed travel with 
different yet constant subsonic speeds, indicating subsonic 
radial oscillations; (2) the first shock near the core is an ac- 
cretion shock, as the reduced radial speeds of the pre-shock 
and post-shock are both negative for inflows. This accretion 
shock expands at a constant subsonic speed. For the same 
x*(2) and x t , the velocity of the first shock in Class I solu- 
tion is a bit larger than that in Class II solution, and for the 
same class of shock solutions, the velocity of the first shock 
becomes larger as the a;* (2) becomes smaller; and (3) the 
location of the second shock x s (2) can be larger or smaller 
than 1, such that the velocity of this shock can be supersonic 
or subsonic. The upstream solution for this shock evolves 
from the accretion condition to breeze to wind. For different 
astronomical system, we can construct different shock solu- 
tions (e.g., an envelope contraction with core collapse or an 
envelope expansion with core collapse). 
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Table 4. Class I and II Shock Breeze Solutions with Different r Values 



Parameters 


T 


A 






Vd 


ot-d 


Vu 




x,(l) = 0.23 


1 


1.86 


1.37 


1.37 


0.529 


1.54 


0.181 


1.09 




1.3 


2.66 


1.34 


1.74 


0.515 


1.58 


-0.450 


0.770 




1.5 


3.31 


1.31 


1.62 


0.500 


1.62 


-0.731 


0.732 




1.8 


4.49 


1.27 


2.29 


0.480 


1.69 


-1.12 


0.704 


B = 0.10 


1 


0.23 


2.42 


2.42 


1.76 


0.173 


0.898 


0.0746 




1.3 


0.48 


2.33 


3.03 


1.69 


0.166 


0.575 


0.0564 




1.5 


0.58 


2.29 


3.44 


1.66 


0.164 


0.452 


0.0521 




1.8 


0.78 


2.25 


4.05 


1.62 


0.161 


0.316 


0.0485 



Columns 1 to 9 contain shock solution properties [x»(l) for Class I solutions 
and B for Class II solutions]; the ratio r of the downstream (a^) to upstream 
(a u ) sound speed; the mass parameter A for the upstream solution; the down- 
stream shock location x st i\ the upstream shock location x su ; the downstream 
reduced velocity Vd', the downstream reduced density ay; the upstream re- 
duced velocity v u and the upstream reduced density a u , respectively. 



Table 5. Class II Shock Solutions for r = 1.5 with B = 10" 6 . 







Vd 


Old 








A 




V 


2.0 


3.0 


1.44 


1.58 x 10~ 6 


-0.209 


4.12 X 10" 


7 


8.41 x 10" 


-7 


-0.954 


3.0 


4.5 


2.30 


2.73 x 10~ 6 


1.65 


1.00 x 10" 


6 


4.65 x 10" 


6 


0.998 


3.5 


5.25 


2.76 


3.79 x 10- 6 


2.47 


1.51 x 10- 


6 


9.61 x 10- 


6 


1.85 



Columns 1 to 8 are the downstream shock location x s( j; the upstream shock location 
x su ; the downstream reduced radial speed the downstream reduced density ay; 
the upstream reduced radial speed v u ; the upstream reduced density a u and the mass 
parameter A and the speed parameter V for the upstream solution, respectively. 
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Figure 16. The negative reduced radial speed —v(x) versus x 
for r > 1. The straight dashed line is the sonic critical line. The 
light solid curves form the Class I shock breeze or contraction 
solutions with x» =0.23 and different r > 1, and the heavy solid 
curves form the Class II shock breeze solutions with B = 0.10 for 
different r > 1. Here, we choose r to be 1.0, 1.3, 1.5 and 1.8 to 
get different shock breeze or contraction solutions. 



3.5 Two- Temperature Similarity Shock Solutions 

In the preceding sections, we only consider the case of r = 1 
corresponding to the same thermal temperature for both 
downstream and upstream fluids such that the entire fluid 



Figure 17. The negative reduced radial speed — v(x) versus x. 
The straight dashed line is the sonic critical line. The solid curves 
form the Class II shock solutions with B = 10~ 6 for r = 1.5. We 
choose x s d to be 2.0, 3.0, 3.5 to construct different shock solutions 
matched with asymptotic inflow and winds. 



has the same sound speed a u = ad- For the isothermal con- 
dition, the sound speed can be expressed as 



1/2 



(Z+l)fc J V2 



(23) 



where k is the Boltzmann constant, and sound speed a is 
determined by the mean atomic mass p, the ionization state 
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Z and the temperature T together. In some astrophysical 
processes, these factors determining the sound speed should 
evolve from an interior region to an exterior region (e.g., Hll 
regions around luminous massive OB stars); in this sense, 
a fluid may not be treated as an isothermal system. In this 
section, we will discuss a system of two-temperature fluid. 
In such a system, the downstream sound speed ad differs 
from that of the upstream a u . Meanwhile, the downstream 
and upstream fluids can be treated as isothermal fluids re- 
spectively. In most astrophysical systems, the downstream 
sound speed should be larger than that of the upstream. As 
the central objects heat up the gas around them, we presume 
t = ad/a u > 1 in our numerical examples. 

Let us begin with the shock solutions for the case of 
V — at x — > +oo. Here, we take on Class I solutions with 
a;*(l) = 0.23 and Class II solutions with B = 0.1 as exam- 
ples of illustration. By matching the upstream and down- 
stream solutions in the a — v phase diagram at a meeting 
point xf, we obtain the Class I and Class II shock solu- 
tions for different values of r (e.g., r = 10, 1.3, 1.5, 1.8). In 
Fig. 1161 we display both classes of shock solutions for — v(x) 
versus x. Relevant parameters are summar ized in Table H 
For a prescribed boundary condition at x — » + [i.e., mo 
corresponding to x*(l) or the central reduced mass density 
B], as parameter r increases, the downstream shock loca- 
tion x s d representing the Mach number of downstream shock 
becomes smaller; meanwhile, the mass parameter A of the 
upstream solution increases. For A < 2, the upstream so- 
lution gives an outflow breeze; as the value of A increases, 
the breeze becomes weaker. For A > 2, an inward contrac- 
tion appears at x — > +oo; the strength of this contraction 
or inflow becomes larger, as the value of A increases. 

We then allow the velocity parameter V 7^ to vary. 
As an example, we take r = 1.5 to derive the Class II shock 
solutions. By choosing different downstream shock locations 
x S d, we construct and display different shock solutions for 
—v(x) versus x in Fig. 1171 with the reduced central density 
B — 10 -6 . Details of such type of two-temperature shock 
solutions are summarized in Table |5] From this family of 
shock solutions, we note that with increasing value of x s d, 
the upstream solutions gradually change from an accretion 
to a wind, and the shock strength becomes weaker. As x s d is 
increased further, Vd will becomes smaller than v u (e.g., x s 
1 for t — 1.5) which is unphysical. That is, the downstream 
shock location x s< j cannot be arbitrarily large. 



4 ASTROPHYSICAL APPLICATIONS 

Self-similar solutions have been studied to model various as- 
trophysical systems previously (e.g., Larson 1969; Penston 
1969; Shu 1977; Tsai & Hsu 1986; Shu et al. 1987, 2002; Lou 
& Shen 2004; Shen & Lou 2004). There are many instances 
in which these solutions allow one to develop a numerically 
simple description of what one first feels is a very compli- 
cated astrophysical flow. As results of instabilities and su- 
personic flows of fluid, accretion and expansion shocks are 
expected to emerge in various cosmological and astrophys- 
ical settings (e.g., supernova explosions; dynamical Hll re- 
gions surrounding luminous massive OB stars; the dynamic 
evolution from AGBs to pPNe; accretion shocks around BHs 
and in quasars; possibly, the creation of initial fireballs in a 



class of GRBs). In this paper, we mainly focus on various 
theoretical possibilities and produce different types of simi- 
larity shock solutions, including the solutions for a static SIS 
envelope, for asymptotic breezes or winds or inflows, for twin 
shocks, for two-temperature fluid and for the global feature 
of EECC solutions as well as various other possible combi- 
nations. Here, we discuss and outline potential astrophysical 
applications for these similarity shock solutions. 

1. The dynamic process connecting the AGB phase and 
the pPN phase (e.g., Lou & Shen 2004). The stellar evolu- 
tionary stage between the end of the AGB phase and the PN 
phase has long been a key missing link in our physical un- 
derstanding for the evolution of a single star (e.g., Balick & 
Frank 2002). A swelling and expanding giant or supergiant 
star with a massive yet slow wind eventually end up with 
a more or less detached system containing a small compact 
and extremely hot white dwarf at its center and a nebula 
of various morphologies with signatures of shocks indicat- 
ing interactions of hot faster winds (~ 10 3 km s" 1 ) with the 
massive slow wind (~ 10 — 20 km s _1 ). In the gross picture, 
a dynamic process characterized by self-similar EECC shock 
solutions may be highly relevant: expansion and outflow of 
gas materials surrounding a proto-white dwarf in a pPNe, 
meanwhile infall and collapse in the central core region to 
produce a compact and hot white dwarf. 

Sufficiently far away from initial and boundary con- 
ditions, a largely isolated stellar system may gradually 
evolve into a dynamic phase in a self-similar manner. Avail- 
able observations suggest the following constraints: (1) the 
timescale of this collapse-expansion stage is estimated to be 
~ 10 3 yrs; (2) the mass of the central white dwarf should 
be less than the Chandrasekhar limit of 1.39 Mq, otherwise 
nova- like or even Type la supernova-like activities might oc- 
cur sporadically to produce anomalous abundances in PNe; 
(3) the massive shell envelope expands with a typical speed 
of ~ 10 — 20 km s _1 which is comparable to the sound speed 
in the gas; (4) there are numerous indications that PNe form 
involving a much slower wind from the progenitor giant star 
overtaken by a subsequent faster hot wind presumably gen- 
erated when a white dwarf emerges, such that shocks are 
inevitable (Kwok 1978, 1982, 1993; Balick & Frank 2002). 

Adopting highly idealized Class I similarity shock solu- 
tions constructed in this theoretical analysis, we attempt to 
provide the following rough estimates in reference to the the- 
oretical central mass accretion rate given by mod 3 /G where 
a is the isothermal sound speed and G is the gravitational 
constant. Here, we take the sound speed a as ~ 20 km s 
with a mass accretion rate of the order of ~ 2mo Mq /yr. 
For a timescale of ~ 10 3 yr, the total accumulated mass onto 
the central white dwarf would be ~ 2 x 10 3 mo Mq . As the 
mass limit of a white dwarf 1.39 Mq and considering the 
possible range of sound speeds, the parameter mo should 
fall in the rough range of ~ 10 -3 — 10 -4 . In our analysis, 
mo associate with the Class I shock solutions with only one 
stagnation point appears to fit in this range. For example, 
we may take a Class I twin shock solution with x*(l) = 0.8 
with a corresponding mo = 6.20 x 10 -4 ; this would give a 
mass of ~ 1.2 Mq for the eventual white dwarf. By choosing 
the second shock location at x s (2) — 1.15, we can derive the 
expansion speed of the second shock to be ~ 20 km s _1 and 
the wind speed of ~6 km s _1 , consistent with observations 
of wind speed in the range of ~ 5 — 10 km s _1 . We note that 
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there exists an accretion shock in the interior region [i.e., 
x s (l)], which expands at a constant speed of ~ 13 km s _1 ; 
this expansive shell form part of interior structures of a PN. 

2. Cloud collapses in the star formation process. Shu 
et al. (1987) reviewed star formation processes in molecular 
clouds. They suggested to divide the star formation process 
into four stages. In this scheme, the self-similar evolution 
may be applicable in relatively early stages referred to as 
the pre-stellar stage of star formation defined as the dynamic 
phase in which a gravitationally bound core has formed in a 
molecular cloud and evolves towards higher degrees of cen- 
tral condensation, but no central hydrostatic protostellar ob- 
ject exists yet (e.g., Andre et al. 1999). As idealization and 
simplification, a molecular cloud may be treated as spheri- 
cally symmetric and isothermal during that stage. 

With this, the self-similar shock solutions to be used in 
star formation processes can be classified into two types: 

(i) The solutions with an asymptotic radial speed u — > 
at r — * +oo. This type of solutions has been discussed by 
Shu (1977) and Tsai & Hsu (1995) for a static SIS enve- 
lope with V — and A = 2 in In § 3.3, we derive 
various breeze solutions with velocity and density profiles 
and different values of mo. For the central mass accretion 
rate M = m a 3 /G, the value of m is 0.975 for the EWCS, 
while in our shock breeze solutions, the range of — 0.975 
for mo is much wider. This is an important theoretical lee- 
way to account for various central mass accretion rates in 
protostars. For a typical temperature of 10 K in molecular 
clouds, we take a sound speed a = 0.2 km s _1 and a cen- 
tral mass accretion rate of M ~ 2 x 10 _6 mo Mq /yr. As 
the timescale of this stage is about 10 6 — 10 7 yrs (e.g., Jes- 
sop & Ward- Thompson 2001), the formation of a protostar 
with a mass of ~ 1 Mq would imply a parameter mo > 0.1. 
Meanwhile as the gravitational energy is released during the 
core collapse, one would expect a point source of infrared 
luminosity Lir ~ a 6 tmo 2 / (GR t ) emanating from the dense 
gas cloud surrounding the core, where i?* is a protostellar 
radius taken to be ~ 3 Rq (e.g., Stahler 1988). Note that 
a factor of m§ is involved here. For mo ~ 1 as in the case 
of an EWCS of Shu (1977), this ml factor is fairly close 
to 1. For smaller values of mo, this m§ factor may reduce 
the Lir, significantly. For the dense core L1544 in the Tau- 
rus molecular complex with a lifetime of ~ 0.3 — 0.5 Myr 
(e.g., Tafalla et al. 1998) and at the distance of L1544, the 
sensitivity of the Infrared Astronomy Satellite (IRAS) was 
about 0.1 Lq . As the LI 544 cloud system could not be de- 
tected by the IRAS, the infrared luminosity Lir should be 
less than 0.1 Lq . From this, we may infer the upper limit 
of mo to be ~ 0.1 for the core of the L1544 cloud system. 
Among other factors and considerations as well as possi- 
bilities, our shock solutions here properly adapted should 
be applicable to L1544 system. We note that Myers (2005) 
recently presented the models of free-fall gravitationl col- 
lapse in the equilibrium layers, cylinders and Bonnor-Ebert 
spheres, and the spherical models can match with the ob- 
served inward velocity very well. Because of the very little 
accretion rate of the collapsing Bonner-Ebert sphere in the 
early stage of the star formation, the model of Bonnor-Ebert 
spheres may be an alternative way of solving the luminosity 
puzzle of LI 544 cloud system involving centrally condensed 
collapse of a starless core. 

(ii) The shock solutions with a constant asymptotic ra- 



dial speed uasr^ +oo. Recent observations have revealed 
that at sufficiently large distances away from the core, the 
envelopes of some molecular clouds may expand with a speed 
of ~ 0.25 km s~\ hinting at the characteristic feature of the 
EECC shock solutions; these theoretical possibilities were 
first discussed by Shen & Lou (2004). 

3. The ionization and expansion of Hll regions. As mas- 
sive OB stars turn on their central nuclear reactions in 
modecular clouds, the interstellar medium of surrounding 
neutral hydrogen gas strongly absorbs the ultraviolet radi- 
ation from OB stars and becomes ionized to form Hll re- 
gions (e.g. Stromgren 1939; Osterbrock 1989 and extensive 
references therein). As the ionization front sweeps through 
the cloud, flows driven by pressure gradients develop be- 
tween Hll and Hi regions as well as within Hll regions, and 
shocks would conceivably emerge in these processes. Phys- 
ically similar processes but on much larger scales are also 
expected to occur in distant quasars and starburst galaxies 
(e.g., redshifts z > 6) where surrounding neutral hydrogen 
gas clouds are irradiated by ultraviolet photons from central 
accretion or starburst activities and are ionized with subse- 
quent dynamical consequences. For distant quasars, features 
in the spectral range between the Gunn-Peterson absorption 
trough (Gunn & Peterson 1965; Scheuer 1965) blueward of 
the Lya emission line (e.g., Fan et al. 2001) and the Lya 
emission line should contain valuable diagnostic informa- 
tion of the underlying dynamics of Hll regions. We shall 
further pursue this line of research by combining radiative 
transfer processes with various shock flow models in sep- 
arate papers. The expansion of the Hll regions has been 
studied for decades (e.g., Newman & Axford 1968; Mathews 
& O'Dell 1969; Tenorio-Tagle 1979; Franco, Tenorio-Tagle 
& Bodenheimer 1990; Shu et al. 2002; Shen & Lou 2004). 
The last two analyses invoked the self-similar shock solu- 
tions to model flows and shocks in Hll regions. Because of 
the heating from the central object and the ionization of 
the gas medium, the sound speed in the inner and outer re- 
gions should be different in general. In our relatively simple 
model framework, a two-temperature fluid separated by a 
shock may grab certain gross thermal features better. For 
a given central reduced mass density, by choosing a sound 
speed ratio of a sl i and a su (e.g., r = 1.5) and the shock lo- 
cation x s d (x S u = TX 3 d), it is possible to construct different 
shock solutions to match with different asymptotic solutions 
(e.g., outflow, inflow and static envelope) at x — » +oo. 

We note that when the isothermal approximation is re- 
placed by the polytropic approximation (Suto & Silk 1988; 
Lou & Gao 2005 in preparation; Wang & Lou 2005 in prepa- 
ration), it would be more natural to deal with temperature 
variations in flows and across shocks. 

4. Accretion shocks around supermassive black holes 
such as quasars. Gas materials around quasars fall towards 
the massive centre following the immense pull of gravity. As 
the infalling gas with a speed of ~ 500 km s" 1 impacts on 
the core gas materials, a strong accretion shock emerges. 
Downstream of such a shock, neutral hydrogen gas is ex- 
pected to be fully ionized. Because of absorptions by the 
infalling pre-shock gas, a sharp radiative flux drop is ob- 
served in quasar spectra indicating the existence of accre- 
tion shocks. For the conceptual simplicity, the geometry of 
such an accretion shock is approximated as spherically sym- 
metric. The temperature of surrounding shock heated gas is 
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higher than ~ 10 K with a sound speed of about several 
10 2 km s _1 . Part of the gas is expected to collapse onto the 
galactic disk eventually. On larger scales, we may apply our 
similarity shock solutions to model accretion shocks for such 
a quasar system. In our scenario, the upstream flow velocity 
is ~ 0.5 sound speed corresponding to several 10 2 km s _1 
and the shock front expands at a constant speed of several 
10 2 km s _1 . As the lifetime of a quasar is about 10 8 yrs, 
the radius of such a shock sphere is about 0.1 Mpc and the 
central mass accretion rate can reach 10 2 Mq yr _1 . This 
accumulation of materials can assemble the host galaxy of a 
quasar in a timescale of several 10 8 yrs. These estimates are 
well consistent with the results of Barkana & Loeb (2003). 
Generally speaking, in the reservoir of our shock solutions, 
we could construct other shock solutions that may match 
with the quasar spectra better. By adjusting locations of the 
second shock, it is possible to match with various asymptotic 
envelope solutions. 

5. Supernova (SN) explosions and evolution of SN rem- 
nants. During the late evolution stage of a massive star, 
the photodisintegration and the neutronisation of the iron 
core will reduce the pressure of the inner stellar core, trig- 
gering off a rapid core collapse. When the inner core stops 
collapsing and rebounds the infalling materials, shocks will 
form and propagate outwards through the inner core, stellar 
interior and envelope. As the SN envelope expands during 
the explosion of a SN, it will interact with the circumstellar 
medium (CSM) and drive another shock through the pow- 
erful stellar wind. We note that the speed of SN ejecta could 
be ultra-relativistic (e.g., Blandford & McKee 1976; Ostriker 
& McKee 1995; Bisnovatyi-Kogan & Silich 1995) and a gas 
may be modelled as polytropic (e.g., Goldreich & Weber 
1980). In this sense, our isothermal shock solutions may not 
be directly applicable. However, the relevant concepts and a 
variety of possibilities shown in this paper are expected to be 
work for a polytropic gas as well. For example, a shock can 
not only match with a static envelope but can also match 
with a wind, as the CSM is contributed to stellar winds in 
the late evolution stage of a massive star. 

By discussions above, we show that similarity shock so- 
lutions are potentially applicable to various astrophysical 
systems, ranging from stellar to galactic scales. There exist 
collapses, outflows or both collapses and outflows with one 
or more outward moving shocks in these systems. 



5 SUMMARY AND CONCLUSION 

In this paper, we explore the self-similar isothermal shock 
solutions, in reference to the earlier work of Tsai & Hsu 
(1995), Shu et al. (2002) and Shen & Lou (2004). 

First, we have significantly expanded the Class I and 
Class II shock solutions of Tsai & Hsu (1995) to two families 
of infinitely many discrete shock solutions matched with a 
static SIS envelope by systematically searching for intersec- 
tions in the a — v phase diagram. These shock solutions have 
subsonic radial oscillations (Shen & Lou 2004). As the value 
of reduced central mass mo decreases, the number of stag- 
nation points in the corresponding v(x) solution increases. 
For a given mo, the location x, of the solution crossing the 
sonic critical line is fixed. The shock in the solution with two 
stagnation points is an accretion shock. It is also possible to 



construct various similarity shock solutions with different 
values mo to match with asymptotic breeze solutions. 

Secondly, we further obtained shock solutions with twin 
shocks matched with a static SIS envelope. This type of 
shock solutions contains an inner accretion shock and an 
outer expansion shock. Both shocks travel outward with con- 
stant yet different speeds. The travel speed of the accretion 
shock is about 0.6 times the sound speed and that of the 
expansion shock is about the sound speed. By adjusting the 
outer expansion shock location, twin shock solutions can 
match with various asymptotic solutions [e.g., inflow (ac- 
cretion) or outflow (wind)]. We apply this class of shock 
solutions to bridge the AGB phase and pPN phase. By the 
property of EECC shock solutions, we suggest to model ac- 
cretion processes with shocks around black holes. 

Finally, by allowing different yet constant temperatures 
on two sides of a shock (i.e., r > 1 with the downstream tem- 
perature higher than the upstream temperature) , we derived 
shock solutions in a two-temperature fluid to match with 
various asymptotic solutions at x — > +oo. These shock solu- 
tions may describe cloud collapse leading to star formation 
and expansion of Hll regions surrounding massive OB stars. 

In our model analysis, we pursue a simple isothermal 
scenario. This is a special case of the more general poly- 
tropic treatment (Cheng 1978; Goldreich & Weber 1980; 
Yahil 1983; Bouquet et al. 1985; Suto & Silk 1988; Maeda 
ct al. 2002; Harada et al. 2003; Fatuzzo et al. 2004; Wang 
& Lou 2005 in preparation). Shock solutions have also been 
derived in the polytropic case (Lou & Gao 2005, in prepara- 
tion) . Polytropic shock solutions should have an even wider 
range of astrophysical applications. 
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Figure Al. The phase diagram of v versus a at the meeting 
point xp = 0.5. Each asterisk symbol denotes an integration back 
from the shock location x s with its value shown explicitly; each 
solid dot denotes an integration from x — * + for an inner LP 
solution [equation llll l with the parameter B marked explicitly; 
each open circle denotes an integration from x* with its value 
shown explicitly and this integration starts from a type 2 eigen- 
solution across the sonic critical line x — v = 1. For Class I shock 
solution, we determine the location x, = 0.0544 where the flow 
solution crosses the sonic critical line analytically and the shock 
location x s = 1.26 from the intersection point of the solid and 
dotted curves in the v — a phase diagram. For Class II shock so- 
lution, we deterimine the reduced core density B = 7.9 and the 
shock location at x s = 1.34 from the intersection point of the 
dashed and dotted curves. 



APPENDIX A: SOLUTION MATCHING 

We determine the parameters x„, x s and B by continuously 
varying relevant parameters to search for intersections of the 
two curves v and a in the phase diagram of v versus a where 
the meeting point is taken to be xf — 0.5. We integrate non- 
linear ODEs J7J and (|HJ from s» (< if) by using the type 
2 eigensolutions across the sonic critical line as specified by 
equation I13H to xf = 0.5 for Class I solutions or indepen- 
dently, from x — > + by using the LP-type solution specified 
by llll l as the starting condition to xf = 0.5 for Class II 
solutions. In this manner, we obtain a pair of v and a in the 
so-called phase diagram. By gradually varying parameter x* 
(for Class I) or B (Class II), the values of v and a at the 
meeting point xf will change continuously in the v — a phase 
diagram. In Fig. IA1I we show such phase curves in solid and 
dashed curves for Classes I and II solutions, respectively. 

Meanwhile, we choose different shock locations and ap- 
ply the jump condition (1201 to obtain compatible Vd and 
ay as the starting condition and integrate back from x 3 to 
xf = 0.5 to determine values of v and a at if = 0.5; the 
dotted curve in Fig. IA1I is obtained this way and from the 
intersections between the dotted curve and the solid and 
dashed curves, we determine parameters x,, x 3 and B for 
relevant similarity shock solutions. 

The mass parameter A and shock location x„ for the 
shock solution of Shu et al. (2002) for a given B can be 
determined from the phase diagram of v versus a. Two spe- 
cific examples are shown in Fig. IA2I for B = 1 and B = 4, 
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Figure A2. The phase diagram of v versus a at the meeting 
point xp = 3. Each asterisk symbol indicates an integration from 
the shock location x s with its value shown explicitly for B = 4; 
each open circle symbol indicates an integration from the shock 
location x s with its value shown explicitly for the case of B = 1; 
each solid dot denotes an integration from x — > +oo (practically, 
we take x = 20) to use the analytical asymptotic solutions lO 
with the velocity parameter V = and the mass density param- 
eter A with its value marked explicitly. From the intersections 
of the dashed curve with the dotted curve and the solid curve, 
respectively, we determine the shock location x a and the value of 
mass density parameter A for B = 4 and B = 1, respectively. 



with the meeting point chosen at if = 3. The numerical 
procedure is similar to that described in § 3.1.1. 

This paper has been typeset from a TfrjX/ I^TgX file prepared 
by the author. 



